

% Find the population cutoff
CumPopC = cumsum(popC) ;
xCut = FindZeroInterp( CumPopC - 0.023 , Grids.x ) ; % 0.023
IsInZ = Grids.x <= xCut ;
imaxZ = sum(IsInZ) ;

% adjust for profit-to-GDP ratio in model vs. data
ProfToGDP = AgProfits(vC,zetaC,GridsC) / AgOutput(vC,zetaC,GridsC) ;
scale     = 0.3 / ProfToGDP ;

% compute phi1Z with a smooth S-shape function
Grids.SmoothPolicy = 50 ; % smoothing parameter: 50 
phi1Z = 1 - normcdf( ( Grids.x - xCut ) * Grids.SmoothPolicy ) ;
tau1Z = 0.0004 / scale * AgOutput(vC,zetaC,GridsC) / AgTaxExp(1,1+phi1Z,vC,zetaC,GridsC) ;
Grids.tax = 1 + tau1Z * phi1Z ;

% compute phi1Z derivative with a smooth S-shape function
Grids.DlogPolicy = - tau1Z * Grids.SmoothPolicy * normpdf( ( Grids.x - xCut ) * Grids.SmoothPolicy ) ./ Grids.tax ;          
    % also adjust by scale to reflect that capital use would magnify the
    % impact of the subsidy with a Cobb-Douglas production function